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Convergence of multigrid and defect-correction iterations is comprehensively studied within different incom- 
pressible and compressible inviscid regimes on high-density grids. Good smoothing properties of the defect- 
correction relaxation have been shown using both a modified Fourier analysis and a more general idealized- 
coarse-grid analysis. Single-grid defect correction alone has some slowly converging iterations on grids of 
medium density. The convergence is especially slow for near-sonic flows and for very low compressible Mach 
numbers. Additionally, the fast asymptotic convergence seen on medium density grids deteriorates on high- 
density grids. Certain downstream-boundary modes are very slowly damped on high-density grids. Multigrid 
scheme accelerates convergence of the slow defect-correction iterations to the extent determined by the coarse- 
grid correction. The two-level asymptotic convergence rates are stable and significantly below one in most of 
the regions but slow convergence is noted for near-sonic and very low-Mach compressible flows. Multigrid 
solver has been applied to the NACA 0012 airfoil and to different flow regimes, such as near-tangency and 
stagnation. Certain convergence difficulties have been encountered within stagnation regions. Nonetheless, for 
the airfoil flow, with a sharp trailing-edge, residuals were fast converging for a subcritical flow on a sequence 
of grids. For supercritical flow, residuals converged slower on some intermediate grids than on the finest grid 
or the two coarsest grids. 


I. Introduction 

Defect correction (DC) is currently a cornerstone approach for solving the Euler and Navier-Stokes equations. 
Second-order finite-volume discretizations (FVD) require large-stencil linearizations, making direct iterations expen- 
sive. Also, linearizations of inviscid discretizations beyond first-order are highly non-positive and difficult to relax. 
On the other hand, upwind-biased first-order equations are more diagonally dominant and can be relaxed (solved) with 
conventional approaches. Thus, DC is widely used for second-order solutions, 1 * either directly by solving a series of 
first-order equations with modified residuals or indirectly by using the first-order operator to relax or precondition the 
second-order equations. The concept is also being applied in p-multigrid methods to solve higher-order discretizations. 

Usually, DC is cited as being slow to converge the second-order residuals but fast to converge quantities of en- 
gineering interest, such as lift and drag. 1-3 On the other hand, DC has been used to solve large-scale turbulent ap- 
plication problems for many years 4-6 and relatively fast asymptotic convergence of residuals has been observed in 
many instances. A hierarchical full-approximation scheme (FAS) multigrid method 6-7 with a DC-based relaxation 
scheme, herein referred to as MG-DC, was previously developed and applied in two dimensions, demonstrating fast 
convergence of residuals for airfoils at compressible and incompressible conditions. 

Analysis of DC convergence for two-dimensional (2D) convection has been previously performed in a semi- 
discrete setting 6- 8 in which boundary conditions in one direction are taken into account. A two-level multigrid analysis 6 
showed that although the number of cycles to attain convergence was dependent on the mesh density, the dependence 
was reasonably small and fast asymptotic convergence was eventually attained. A more detailed study of DC alone 8 
showed that an asymptotic convergence of about 0.5 per DC iteration is observed in computations. Slow convergent 
DC iterations may be encountered for nonaligned flows before attaining the asymptotic rate; the number of slow itera- 
tions slightly grows on finer grids as /i - l/3 , where h is a characteristic mesh size. This h dependence can be observed 
for three-dimensional flows as well. 

With the current trend of performing complex computations on increasingly larger scales, it is critically important 
to (re)evaluate performance of traditional algorithms on grids of high density. Analysis of convergence on such grids 
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has been conducted in this paper. Some surprising results have been obtained regarding the DC asymptotic rate. 
Specifically, the asymptotic convergence on typical computational grids is significantly different from the asymptotic 
convergence on high-density grids. The asymptotic rates are essentially invariant for several refined grids of medium 
density, but the convergence rates slow significantly with progressive grid refinement. The asymptotic slowdown on 
high-density grids was found first for the Euler system of equations, but was found to occur even for the convection 
equation alone. The previous results for asymptotic convergence of DC iterations are revisited in the light of these 
new findings. 

The purpose of this paper is to analyze convergence of iterative solvers for inviscid flows, ranging from incom- 
pressible to supersonic Mach numbers, to complement the methodology developed previously for diffusion. 9, 10 The 
convergence of the MG-DC algorithm is comprehensively studied within different incompressible and compressible 
regimes on structured grids of progressively high density. The approach is to first assess the convergence away from 
any boundaries and discontinuities that may exist and this assessment can be performed using the framework of a 
small-perturbation (SP) flow. With acceptable and quantified performance within this regime, a solid foundation is 
established for assessing convergence for the general 2D inviscid flow. The entire flow field around an airfoil, for 
instance, has at least six distinct regions (regimes): (1) flow away from the boundaries and discontinuities; (2) flow 
near tangency boundaries away from stagnation; (3) flow within the leading-edge (LE) stagnation flow; (4) flow within 
the trailing-edge (TE) stagnation flow; (5) flow near discontinuities , e.g., shocks; and (6) flow near the far boundary. 
Each of these flow regimes may introduce difficulties in the multigrid and each should be studied individually, both 
analytically and computationally. 

Several analysis tools are used to characterize performance of the MG-DC scheme. For SP flows, a constant- 
coefficient approximation is analyzed with the local mode Fourier (LMF) analysis and a semi-discrete (SD) analysis. 
General quantitative analysis tools 10, 11 idealized coarse-grid (ICG) and idealized relaxation (IR), are applied in actual 
flow computations for assessing multigrid relaxation and coarse-grid correction. The analytical results, confirmed with 
actual computations, indicate that asymptotic MG-DC convergence rates are stable and well separated from one and 
are limited on high-density grids by the quality of the coarse-grid correction. The convergence of MG-DC iterations 
is significantly better than convergence observed in DC iterations alone because multigrid accelerates convergence of 
slow DC iterations, especially for near-sonic flows and low-Mach compressible flows. 

The material in the paper is presented in the following order. Components of the multigrid and defect correction 
scheme are presented in Section II. Analysis tools are introduced in Section III. Section IV describes the first-order 
solver that serves as a driver of the DC iterations. An analysis of DC and MG-DC iterations for SP flows is presented 
in Section V. Numerical tests and IR/ICG analysis of flows in other regimes are presented in Section VI. The 
results are discussed in Section VII. Details of the analysis methods used in this paper are provided in Appendices A 
and B. Asymptotic convergence rates of constant coefficient convection equation on high-density grids are discussed 
in Appendix C. For reference. Table 1 includes all acronyms used in the paper. 

II. Components of MG-DC solver 

The conservation form of the 2D steady inviscid flow equations is given as 

R(Q) = 0. (1) 

Here, the conserved variables for compressible flows are Q = (pu, pv, pw, p, pE ) T , representing the momentum 
vector, density, and total energy per unit volume, and R(Q) is a spatial divergence of convective fluxes 

R(Q) = <9 X F(Q) + d y G(Q), (2) 
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The primitive flow variables are velocity, pressure, and density, q = (it, v, p. p) T . Eq. (1) is discretized with a second- 
order, cell-centered, upwind-biased FVD scheme that employs an approximate Riemann solver to compute fluxes at the 
control volume faces. The baseline Riemann solver is the flux-difference-splitting (FDS) scheme 12 but other schemes 
are also considered, including the low-dissipation flux-splitting (LDFS) 13, 14 and flux-vector-splitting (FVS). 15 Either 
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Acronym 

Description 

Alternating 
Line-Colored (ALC) 

A relaxation method 

Courant-Friedrichs-Lewy 

An iterative parameter characterizing the ratio of 

(CFL) number 

(pseudo) time increment to mesh spacing 

Correction Scheme (CS) 

A MG scheme that uses linear approximations on coarse grids 

damped Alternating 
Line-Jacobi (dALJ) 

A relaxation method 

Defect Correction (DC) 

Used for single-grid iterations and as relaxation in multigrid 

Flux Difference 
Splitting (FDS) 

A less-dissipative approximate Riemann solver 

Flux Vector 
Splitting (FVS) 

A more-dissipative approximate Riemann solver 

Full-Approximation 
Scheme (FAS) 

A MG scheme that uses non-linear approximations on coarse grids 

Full-Multigrid 

A MG scheme that uses coarser-grid solutions 

(FMG) 

to form finer-grid initial approximations 

Finite-Volume 
Discretization (FVD) 

The discretization approach used in this paper 

Idealized 

Coarse Grid (ICG) 

General quantitative method for analysis of multigrid relaxation 

Idealized 
Relaxation (IR) 

General quantitative method for analysis of coarse-grid correction 

Leading Edge (LE) 

Designate the leading-edge stagnation area 

Low-Dissipation 
Flux-Splitting (LDFS) 

A more-dissipative approximate Riemann solver 

Local-Mode 

A constant-coefficient analysis for interior of the domain. 

Fourier (LMF) 

assumes periodicity in all directions 

Multigrid (MG) 

A hierarchical computational method 

MG-DC 

Multigrid method studied in this paper 
that uses defect-correction based relaxation 

Semi-discrete (SD) 

A constant-coefficient analysis taking boundary conditions into account, 
assumes periodicity in the directions tangential to the boundary 

Small Perturbation (SP) 

Computational model that assumes small deviation 
from a known (e.g., free-stream) solution 

Trailing Edge (TE) 

Designate the trailing-edge stagnation area 


Table 1. Acronyms used in this paper. 


of these latter schemes are generally known to be more dissipative than the FDS scheme. The discrete approximations 
to derivatives correspond to the Fromm discretization for structured grids used herein. 

The same approach is used for incompressible flows with small variations. The variables are Q = (it, v, p) T and 
the fluxes are defined as in Eq. (2), except the density is constant and the fourth (energy) equation is dropped. The 
incompressible version of the FDS scheme 4 is used. 

In a DC relaxation, a correction, <5Q^, to the approximate solution, Q h , is computed from the driver equation 

DJQ' 1 = — R(Q fc ), (3) 

where D is the Jacobian of the first-order upwind discretization, and R is the discretized residual Eq. (1). For DC 
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relaxation within an outer FAS multigrid cycle, a correction scheme (CS) multigrid is applied to determine SQ h . One 
CS cycle generally reduces the residual of Eq. (3) by an order of magnitude (see Section IV). For DC iterations, Eq. (3) 
is solved to high precision. Where practical, e.g., for the SD analysis or a scalar convection equation, Eq. (3) is solved 
precisely, otherwise multiple multigrid cycles are used. 

After DC relaxation, the solution of the target FVD scheme is updated as 

Q h = Q h +6Q h . (4) 

For the MG-DC solver, FAS multigrid is used to accelerate convergence. An F AS{y \ , v 2 ) multigrid cycle starts on the 
target finest grid, performs v\ relaxations on the current grid, restricts solutions and residuals to the coarser grid, solves 
the coarse-grid problem recursively, prolongs the coarse-grid correction, and completes with additional z/ 2 relaxations. 
Each coarse grid is obtained by full coarsening from the finer grid. The same FVD scheme is used on all grids and W 
cycles are used. For SD computations, the restriction operator is full weighting, and the prolongation operator is the 
normalized transposition to the restriction. For fully discrete computations, the restriction operator is the conservative 
residual restriction and prolongation corresponds to linear interpolation. Full multigrid (FMG) requires a high-order 
prolongation for full efficiency. In the current FMG solver, the FMG prolongation is the same as within the FAS cycle. 

III. Analysis tools 

In recent years, a number of powerful methods have been developed to analyze convergence of iterative solvers. For 
problems well described in terms of small perturbations, e.g., SP flows, analysis of a constant coefficient approximation 
allows one to estimate various convergence characteristics, such as stability, asymptotic and maximum convergence 
rate, number of slow iterations, etc. For more general problems, windowing and downscaling techniques 16 can be used 
to analyze accuracy and grid convergence of discrete solutions. Quantitative analysis methods, IR(^i, z/ 2 ) and ICOf iq , 
z/ 2 ), l0 ‘ 1 have proved to be invaluable for assessing components of multigrid solvers for general problems. Here, 
parameters and z/ 2 denote the respective numbers of fine-grid relaxations performed before and after coarse-grid 
correction. 

Ill A. Analysis of constant coefficient equations on regular grids 

A constant coefficient linearization to the FVD schemes used here on Cartesian grids is given by 

A+d~w h + A~d+w h + B +d~w h + B~d+w h = 0, (5) 

where w h is a discrete solution vector. For compressible flow, the variables are taken following Mulder 17 as w h = 
( Su,5v,6p/(pc),SS ) T , c is the speed of sound, and S = log (p/p 1 ) is the specific entropy. For incompressible 
flow, w h = ( Su , Sv, Sp) T . The operators, 9“ and d~ are upwind discretizations of derivatives, and 9+ and 9+ are 
downwind discretizations of derivatives. 

Different linearizations are associated with each splitting scheme. For the baseline FDS scheme, the linearizations 
are eigenvalue splittings of the Jacobian matrices associated with non-conservative formulations, 

A = A + +A , B = B+ + B , 

where 
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The speed of sound is taken as c = 1 and the velocities are defined as u = M cos(a), v = M sin(a), where M is 
Mach number and a is the angle of attack. 

For subsonic regimes. 



/ 

{u + c)/2 

0 

(u + c)/2 

0 \ 


^ (u — c)/2 

0 

— (u — c)/ 2 

° \ 

A+ = 


0 

u 

0 

0 

, A- = 

0 

0 

0 

0 



(u + c)/2 

0 

(u + c)/2 

0 


— (u — c)/2 

0 

(u — c)/2 

0 


l 

0 

0 

0 

u ) 


l 0 

0 

0 

0 / 


4 of 17 


American Institute of Aeronautics and Astronautics 



B 


( V 

0 

0 

° ) 


( o 

0 

0 

° ) 

0 

( v + c)/2 

(ti + c)/2 

0 

B = 

0 

{v - c)/2 

-( v - c )/ 2 

0 

0 

(v + c)/2 

(v + c)/2 

0 

5 

0 

-0 - c)/2 

(v - c)/2 

0 

V o 

0 

0 

v J 


l 0 

0 

0 

0 / 


For supersonic regimes, 

A + = A; A - = 0; B + = B; B~ = 0. 

The discretization is defined as 

A + d~w h + A~d+w h + B + d~w h + B ~d+w h = 0. 

In the LMF analysis, the iterations are considered on a periodic domain for discrete Fourier components w h = 
exp (i(8 x i x + Oyiy)), where i x and i y are integer grid indexes. The Fourier frequencies are normalized: \9 X \ < 
7 r, 0y < 7 r. The outcome of the Fourier analysis is an iteration symbol, which is a 4 x 4 matrix with complex 
coefficients parametrized by the Fourier frequencies. The specific grid size is reflected through the range of Fourier 
frequencies realizable on the given periodic grid. For elliptic equations, the maximum spectral radius of the LMF- 
symbol matrix, taken over all realizable frequencies, is an accurate indicator of the asymptotic convergence rates. For 
non-elliptic equations, the spectral radius of the LMF symbol is not a sharp estimate for the asymptotic convergence 
rate on grids of moderate sizes because the LMF analysis accounts only for local error damping, but does not account 
for boundary effects and error propagation along the characteristics. Note, however, that LMF analysis provides a 
useful stability test. A larger-than-one LMF spectral radius is an indication of unstable iterations. 

For multigrid computations, the relaxation smoothing rate is an important characteristic. The smoothing rate is 
estimated as the maximum spectral radius of the LMF relaxation symbol, where the maximum is taken over high 
frequency modes. Typically, high-frequency modes are defined as the modes with max(\0 x \, \8 y \) > all other 
modes are considered smooth. A more general approach is to define the high-frequency modes as the modes that have 
relatively large contributions to the residual. 19 An implication of this definition for non-elliptic problems is that the 
typical set of high-frequencies is reduced: the modes that are smooth in the characteristic directions are excluded, 
even if their Cartesian frequencies are high. For illustration, for the convection flow at 45° discretized on a uniform 
Cartesian grid, the mode exp (i(9 x i x + 9 y i v )) with 9 X ss 7r and 9 y « — 7r is not a high-frequency mode because 
9 X + 9 y ~ 0, and the mode is constant along the characteristic direction. With this modification for non-elliptic 
problems, the LMF analysis predictions of the smoothing rate are reasonably accurate. A more detailed description of 
the modified LMF smoothing analysis is provided in Appendix A. 

The SD analysis is a good predictor of the asymptotic convergence for non-elliptic problems. The SD analysis 
assumes solutions in the form w h = exp(i9 y i y )W h (i x ) , i.e., the solution is a product of a Fourier component in the 
//-direction and a discrete function, W h (i x ), representing solution variations in the a’-direction. The SD analysis is 
accounting for boundary effects and error propagation along the characteristics. For each //-directional Fourier fre- 
quency, the asymptotic rate is estimated as the spectral radius of the SD iteration matrix, which has a size proportional 
to the number of degrees of freedom in the .r-dircction. Another useful feature of the SD analysis is the capability 
to identify slow convergent iterations, characterize the error components causing the slow convergence, and explain 
the mechanism of transition from the slow intermediate convergence to good asymptotic convergence. A detailed 
description of the SD analysis is provided in Appendix B. 

III.B. General quantitative analysis 

More general, quantitative analysis methods for multigrid solutions are IR and ICG iterations. The iterations are de- 
signed to identify slow relaxation or inefficient coarse-grid correction of a multigrid solver. In these iterations, one 
part of the cycle (coarse-grid correction for IR iterations and relaxation for ICG iterations) is actual and its compli- 
mentary part is replaced with an ideal imitation. The IR and ICG methods can be applied to any formulation with a 
manufactured solution; typically zero solution is used. The initial solution is chosen randomly. In IR iterations, the 
relaxation in the cycle is replaced with an explicit error averaging procedure. In the IR methods used for this paper, 
the error at a node is averaged from all the edge-connected neighbors. ICG cycles use actual relaxation scheme and 
emulate the coarse-grid correction by, first, averaging algebraic errors to the coarse grid and, then, interpolating the 
averaged error back to the fine grid as a correction. The results of this analysis are not single-number estimates; they 
are rather convergence patterns of the iterations that may either confirm or refute expectations indicating what part of 
the actual solver should be improved. 
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The IR and ICG iterations can be directly applied in the most complicated situations including highly variable (or 
nonlinear) coefficients, complex geometries, and unstructured grids. The generality of the analysis makes it a very 
valuable tool for analyzing complicated large-scale computational problems, where no other analysis methods are 
currently available. Properties and specific implementations of IR and ICG methods can be found elsewhere. 10 - 1 1 


IV. First-order Euler solver 



Figure 1. Combined stagnation and tangency grids. 



Figure 2. Number of iteration required to reach machine-zero residuals for the first-order solver. 


Mulder 217 developed efficient 2D multigrid solvers for the first-order upwind discretizations of the inviscid flow 
equations using both full-coarsening and semi-coarsening approaches. He analyzed many relaxation schemes us- 
ing a 2-level LMF analysis and showed that the problem of alignment could be addressed uniformly with damped 
alternating-line-Jacobi (dALJ) relaxation within a full-coarsening framework or with point-implicit relaxation within 
a semi-coarsening framework. In this paper, full-coarsening is used with an alternating-line colored (ALC) relaxation. 
An under-relaxation factor, oj = 0.8, is needed to effectively smooth high-frequency error. 111 The performance of a 
two-color ALC relaxation is similar to performance of dALJ relaxation. 

To illustrate the performance of iterative solvers, computations are performed on a domain around a cylinder. A 
typical grid is the union of the two cylindrical grids shown in Fig. 1, has local near-unity aspect ratios, and spans 
180° of arc sector. Inflow/outflow boundary conditions are applied at all boundaries. The initial solution is a random 
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perturbation of the uniform free stream conditions. 

Fig. 2 compares the number of fine-grid relaxations required to reach the machine-zero residual for single-grid 
ALC relaxations (w = 1.0) and FAS( 2, 1) multigrid W-cycles. 1 One ALC relaxation is counted as two relaxations 
and one W-cycle is counted as six relaxations. Results are shown for the FDS scheme on two grids for a range of 
Mach numbers and for incompressible flow (M = 0). The number of single-grid relaxations approximately doubles 
as the grid is refined by a factor of two in each direction, as expected. The required number of relaxations is highest 
at M ss 1. The number of relaxations is lowest for the higher Mach numbers and, somewhat unexpectedly, for the 
least compressible Mach number of M = 0.01. The number of fine-grid relaxations observed within MG-DC solver 
to reach the same residual tolerance is relatively insensitive to variations of Mach number or grid size. Although not 
shown, the asymptotic convergence per cycle is between 0.2 and 0.4 for all Mach numbers on both grids. 


V. Multigrid for small-perturbation flows 

A previous study showed that, even when the asymptotic convergence rates of DC iterations are fast, a number of 
slow iterations precedes the asymptotic regime. The slow convergence occurs for smooth characteristic error compo- 
nents 18, 19 that are very smooth along the characteristic directions. Such components are removed mainly by accuracy 
propagation from boundaries along the characteristics. Such removal may take many iterations because a less-accurate 
driver propagates cross-characteristic oscillations for only short distances. Eventually, however, the smooth character- 
istic errors are removed and asymptotic convergence is attained. 

In practical computations, the slow DC iterations may be overlooked on relatively coarse grids because the itera- 
tions may arrive to the required solution tolerance before the characteristic components begin to dominate the solution 
error. In order to observe this slowdown, one should carefully choose the initial solution approximation. On finer 
grids, this slowdown is a major factor limiting the solution efficiency. 

Multigrid is expected to accelerate convergence of slow DC iterations. Note that full-coarsening multigrid has 
its own problems with characteristic components. The asymptotic convergence of the characteristic errors in a two 
level cycle can be as slow as 0.75 per cycle 19 because cross-characteristic variations propagate shorter distances on 
coarser grids than on finer grids. The multigrid effects on asymptotic convergence rates are significant only in those 
flow regimes in which the asymptotic convergence of DC iterations is slower than the coarse-grid correction for 
characteristic error components. Such situations occur on fine grids. 



0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 

Mach 



(a) LMF analysis on 128 2 grid 


(b) ICG analysis 


Figure 3. Smoothing rate of DC iterations. 


For a subsequent use in MG-DC cycle, the smoothing rate of DC iterations is estimated with the LMF and ICG 
analysis. Fig. 3 shows the predicted smoothing rates. For all flow conditions (Mach numbers and angles of attack), 
the predicted smoothing rates are excellent and grid independent. The LMF predicts the smoothing rate of between 
0.5 and 0.7, and the ICG predicts the rate between 0.5 and 0.6. The smoothing rates predicted by ICG are slightly 
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better than the rates predicted by the LMF analysis because ICG predicts the reduction of high-frequency errors in a 
multigrid cycle, while the LMF analysis predicts the reduction of high-frequency errors in a relaxation. In general, the 
LMF analysis can be modified to account for the coarse-grid effects on high frequencies. 

The need and benefits of multigrid are illustrated in Fig. 4 by the SD analysis for SP flows. The flow conditions are 
M = 0.3, a = 45°, the y-directional frequency is smooth 0 y = and the initial distribution along the ^-direction 
is random. While the asymptotic convergence for both DC and MG-DC iterations is about the same, around 0.6, the 
slowest convergence rate is significantly slower for DC than for MG-DC iterations. 



Figure 4. SD analysis: convergence of DC and MG-DC iterations on 256 2 grid, M — 0.3, a = 45°, and 9 y = r ‘ \ . 


Fig. 5 shows the asymptotic rates of a two-grid 14(1, 0) with DC relaxation. The rates are computed with the SD 
analysis on two coarse grids. The asymptotic rates of MG-DC iterations are stable and well below unity over most of 
the M — a range. The convergence is slow for near-sonic flows (M « 1) at intermediate angles of attack and for very 
low compressible Mach numbers. 



Mach Mach 


(a) Grid 64 x 64 (b) Grid 32 x 32 

Figure 5. Asymptotic convergence of two-level V(l, 0) cycle computed with the SD analysis. The color bars are scaled by the slowest asymptotic convergence 
observed on the grid. 
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VI. Flows at other regimes 


VI.A. Boundary-tangency and stagnation flows 

A typical grid for computations of flows characterized by boundary-tangency and LE stagnation are shown in Fig. 1 
Stagnation is computed on the most-forward part of the grid and boundary-tangency is computed on the upper-most 
part of the grid; each domain spans 90 deg of arc sector. A compressible-flow manufactured solution is composed 
of the velocities from the exact incompressible cylinder flow along with constant enthalpy and entropy. Medium-size 
grids are considered. The finest grid has 128 cells in both the circumferential and radial directions. Computations are 
shown for FAS(2,1) W-cycles using a maximum of six levels. Inflow/outflow conditions are applied at all boundaries 
away from the cylinder surface. 

Flows characterized by boundary-tangency do not represent difficulties for the MG-DC solver. A typical con- 
vergence history for a series of grids is shown in Fig. 6 for M = 0.3, starting from a random perturbation to the 
exact solution on the left and from FMG interpolations on the right. Starting from random perturbations, the residuals 
converge rapidly in the first cycles, converge more slowly in intermediate cycles, and then asymptotically converge 
faster. Starting from FMG interpolations, the number of cycles needed are considerably smaller and machine-level 
zero residuals are encountered before the faster asymptotic rates are encountered. Although not shown, similar residual 
convergence per cycle is attained with ICG(1,0) and IR(2,1) multigrid cycles. 




(a) Random initial conditions (b) FMG initial conditions 

Figure 6. Residual convergence for boundary-tangency computations with the MG-DC solver; M = 0.3; FDS scheme. 


Within stagnation flows, the Jacobian can differ appreciably from the small perturbation linearization Eq. (5) be- 
cause the contribution for the velocity gradient (e.g., 0(u x )) to the linearization can be comparable with or even 
greater than the contributions from differences in velocity (e.g., 0(u/h)). These terms can subtract from the diagonal 
contributions associated with the momentum equations. For incompressible discretization schemes in which the mo- 
mentum equations can be marched before solving an elliptic equation for the pressure, these velocity-gradient terms 
can cause an error amplification when marching into/from stagnation. 11 Here, we find that similar difficulties arise 
for the MG-DC solver because DC can be unstable. The DC convergence is sensitive to the particular discretization 
schemes used for LE stagnation. For instance, DC does not converge for the FDS scheme but does for the LDFS and 
FVS schemes. Fig. 7 shows convergence of the MG-DC solver for stagnation flow using the FDS scheme (left) and 
the LDFS (right) scheme. An infinite CFL number is used for the LDFS scheme but a CFL of 400 is used for the FDS 
scheme. On the two coarser grids, the MG-DC scheme does not converge for the FDS scheme — a smaller CFL is 
necessary for the scheme to remain stable. On the finer grids, the overall residual convergence of either scheme within 
stagnation is similar to that observed for boundary-tangency computations. 

Although not shown, for TE stagnation, both schemes are unstable without addition of a pseudo time step. Anal- 
ysis of convergence within stagnation leads to a variable-coefficient problem problem that is difficult to analyze using 
LMF analysis. One can devise neighborhoods which provide relevant constant-coefficient approximations to the full 
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linearization, but certain parts of stagnation, such as the stagnation streamline, are inaccessible to a constant-coefficient 
analysis. 11 The stagnation flow analysis was actually a motivating factor for the development of more general quanti- 
tative analysis methods. For the airfoil computations in the next section, we simply use the LDFS scheme. The airfoil 
has a sharp trailing edge which does not seem to cause a problem with this scheme. 



(a) FDS scheme, CFL = 400 (b) LDFS scheme, CFL = oo 

Figure 7. Stagnation computations: residual convergence for the MG-DC solver; finest grid is 128 2 cells, the coarsest grid is 16 2 cells. 


VI.B. Multigrid for airfoil 

Computations for the NACA 0012 airfoil following the Vassberg and Jameson benchmark study 20 are shown here. The 
grid, similar to that used for the study, is generated through a sheared adaptation of a conformal grid around a Karman- 
Trefftz airfoil matching the leading-edge radius and trailing-edge angle of the NACA 0012 airfoil. The grid extends 
150 chords outwards from the airfoil, and has nearly unity-aspect-ratio cells. The second-order accuracy was verified 
in computations with lifting and non-lifting manufactured solutions for the Karman-Trefftz airfoil in incompressible 
flow and in compressible flow at moderate Mach numbers. A compressible-flow manufactured solution was defined 
with the velocities from the exact incompressible Karman-Trefftz solution along with constant enthalpy and entropy. 

Fig. 8 and Fig. 9 shows residual and drag convergence history of FAS(2,1) cycles for the NACA 0012 airfoil at 
subcritical lifting conditions (M = 0.5 and a = 1.25) and supercritical non-lifting conditions (M = 0.8 and a = 0), 
respectively. Six grids were used in the computations. FMG cycles were started on the coarsest grid composed of 16 2 
cells. The finest grid contained 256 cells in the directions around and outward from the airfoil. For the subcritical 
computations, convergence rates per cycle are uniformly fast. Residual convergence per cycle is 0.3 on the finest grid. 
Convergence of drag (and also lift, although not shown) is quite fast, within one FMG cycle. The exact drag is zero, 
reflected in the benchmark level shown as well as the value on the finest grid. The drag is converging with second 
order accuracy although finer grids are necessary to confirm this. 

For the supercritical computation, convergence rates per cycle are quite disparate between grids. The two grids 
before the finest grid in the FMG sequence are converging much slower than the finest grid or the two initial coarser 
grids. No limiter is used in these computations. Drag is again converging within one FMG cycle. The drag is 
converging with second order accuracy to the benchmark level. 

VII. Discussions 

The MG-DC solver used here is similar to the multigrid scheme developed previously. 6 7 The previous scheme 
used alternating-line Jacobi and/or colored relaxations that do not provide sufficient damping of high-frequency errors 
in purely inviscid regions of the flow. Analysis methods were not applied to identify this shortcoming and instead 
other parts of the algorithm were modified to compensate, namely a pseudo-time step limited by a maximum CFL 
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Figure 8. Residual and drag convergence history of FAS(2,1) cycles for NACA 0012 airfoil at subcritical lifting conditions; M = 0.5 and a — 1.25. 


of 0(100) was added to the implicit relaxation operator, relaxation subiterations were performed, and dissipation via 
entropy fixes to all fields was added to the FDS discretization. In the present work, we apply under-relaxation based 
upon optimization of ICG(1,0) cycles, do not add a pseudo-time step except within stagnation, and do not add any 
entropy fixes. 

The convergence of the MG-DC solver has been comprehensively studied within different incompressible and 
compressible inviscid regimes. The properties of the solver away from any boundaries and discontinuities are analyzed 
on high-density grids because this region forms the foundation of the methodology. Within this region, the smoothing 
properties of the scheme have been shown to be bounded away from one using both a modified LMF analysis and 
a more general ICG analysis. DC alone has some slowly converging iterations on grids of medium density. This 
behavior has been shown previously for convection but the convergence for the Euler equations is slower than that 
for pure convection. The convergence is especially slow for near-sonic flows and for very low compressible Mach 
numbers. Additionally, the asymptotic convergence seen on medium-density grids is significantly different from the 
asymptotic convergence on high-density grids. Certain downstream-boundary modes are very slowly damped on high- 
density grids. The FAS multigrid scheme accelerates convergence of the slow DC iterations to the extent determined 
by the coarse-grid correction. The 2-level asymptotic convergence rates are well separated from unity over most of the 
region but slow convergence is noted for near-sonic and very low-Mach compressible flows. 

We have applied the MG-DC solver to the NACA 0012 airfoil and to different flow regimes, such as near-tangency 
and stagnation. The MG-DC solver encounters problems within stagnation regions. The FDS scheme is unstable 
without a time step addition for leading-edge stagnation and all schemes have a problem for smooth trailing-edge 
stagnation. Analysis of the linearization within stagnation predicts difficulties associated with the loss of diagonal 
contributions to the momentum equation linearization within decelerating flow. A pseudo-time step addition can 
provide convergence, although the amount varies from grid to grid. Nonetheless, for the airfoil flow, with a sharp 
trailing-edge, residuals were fast converging for a subcritical flow on a sequence of grids. For supercritical flow, 
residuals converged slower on some intermediate FMG grids than on the finest grid or the two coarsest grids. The 
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Figure 9. Residual and drag convergence history of FAS(2,1) cycles for NACA 0012 airfoil at supercritical non-lifting conditions; M — 0.8 and a = 0. 


cause of the slowdown may be associated with the coarse-grid correction near Mach unity. Also, the lift and drag both 
showed second-order accuracy in grid refinement for subcritical and supercritical conditions. 

A key measure of efficiency for a multigrid method is the number of FMG cycles required to converge algebraic 
errors below the level of discretization errors. Ideally only a single cycle is needed. For both airfoil solutions, algebraic 
errors in lift and drag were well below discretization errors after a single FMG cycle. Another key property for an 
iterative solver is to ensure that the residual can be driven (fast) to the zero level if needed. The MG-DC solver provides 
fast residual convergence. The efficiency of the scheme is limited by the coarse-grid correction. Previous work has 
shown that a modified coarse-grid discretization can substantially improve the correction. The effectiveness of the 
scheme needs to be explored on high-density grids and in the regimes with slower convergence. Local relaxations in 
slow-convergence regions may accelerate convergence even further. 

A. Modified Local Mode Fourier Analysis 

For given Mach number and angle of attack, the respective symbols of the target, T, and driver, D, operators on a 
uniform Cartesian grid with mesh spacing h are defined as 

T(<9 X , 0 y ) = A+ A- (e i0xix + 3 - 5e~ l0xix + e ~ 2i0xix ) 

-A- A- (e~ l6xix + 3 - 5e Wxix + e 2i0xix ) 

_|_B+_L ( e i8yiy _|_ 3 _ 5 e ~Myiy _|_ g-2 i6 y i y ^ ^ ’ 

_B - ^- ( e ~ l0yiy + 3 — 5e i0yiy + e 2i0yiy ) . 

D (e x ,9 v )= A+I(l-e-^)-A-I(l-e^) 

+B+i (1 - e~ i0yiy ) - (1 - e i0yiy ) . 
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The symbol of the DC iteration is a 4 x 4 matrix 


DC(# X , Oy) = I — D -1 T, 


( 8 ) 


where I is the identity matrix. 

The smoothing rate, /i, is estimated as the maximum spectral radius 

p, = max p (DCd) , (9) 

where d is a high-frequency indicator. For a flow with a < 45°, 

d={ 1 ' if max(|6y,|6y) > § & Imod^ + ^ 2?r) - ?r| < f; (1Q) 

) 0, otherwise. 

B. Semi-Discrete Analysis 

The SD analysis considers the solutions in the form of e l6y Wi x , i x = 0, . . . , N x , where N x is the number of 
nodes in the ^-direction. The original multidimensional discrete problem is, thus, translated into a one-dimensional 
problem parametrized by the normalized Fourier frequency, \8 y \ < n. The discrete function Wi x is either a scalar 
solution for the convection equation, or a vector solution for the system of flow equations (5). The analysis takes 
into account specific implementations of boundary conditions and is capable to predict details of solution evolution 
in individual iterations. When zero manufactured solution is used, the round-off error does not affect computations, 
which is critical for the ability to observe asymptotic convergence in computations. SD tests routinely encounter and 
treat residuals as small as 10 -150 . The asymptotic convergence rate can be directly evaluated as the spectral radius of 
the iteration matrix. The analysis is precise for a constant-coefficient formulation with y-periodic boundary conditions. 
A description of the analysis in application to constant-coefficient convection equation is provided in a previous paper. x 

The DC iteration matrix has the form: 


DC = I — D -1 T. (11) 

Here I, T, and D are the identity, target, and driver matrices, respectively. For the convection equation, aw x + 
bw y = /, the matrix T corresponds to the Fromm discretization, with a row composed of the following coefficients: 


T = 


0 — 2 — 

U 47i x 


— 5a 3a ■ j? a 

4 h x ih x 4" -° 2 4/i x 



( 12 ) 


B 2 


b 

Ahy 


(i e i9y + 3 - 5e l9y + e 2l9y ) , 


and the main diagonal coefficient is underlined. D is a driver two-diagonal matrix: 


(13) 


D = 


B\ 



(14) 


Bi = (1 - e~ i£y ) . 

hy 

For the system of equations the corresponding matrices are block diagonal. 
The iteration matrix of a two-level MG-DC V(vi, v 2 ) cycle is 


(15) 


MG = (DC)" 2 CGC(DC)"\ (16) 

CGC = I — PT~ 1 RT. (17) 

Here CGC is the coarse-grid-correction matrix, R and P are restriction and prolongation matrices, respectively, and 
T c is the coarse-grid-operator matrix. The size of the multigrid matrices is twice as large as the size of corresponding 
single-grid matrices because multigrid couples two components corresponding to Fourier frequencies 9 y and 0 y + 7 r. 
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Figure 10. Asymptotic convergence rate of DC iterations for constant-coefficient convection equation computed with the SD analysis. 


C. Scalar convection equation 

The asymptotic convergence of DC iterations for the scalar convection equation is computed with the SD analysis. 
The variation of the asymptotic rate with the grid density and the angle of attack is shown in Fig. 10. The convergence 
plots on grids of moderate size with up to 128 2 degrees of freedom are practically over-plotted. For small angles of 
attack, grids with 256 2 and 512 2 degrees of freedom also show similar rates. On finer grids, however, the convergence 
rates are dramatically different. Slow asymptotic rates are observed for solutions that are exponentially decaying from 
the outflow boundary toward the interior. Fig. 1 1 shows the real and imaginary components of an eigensolution for 
DC iterations on a grid with N x = 2048 and a = 45°; only variation near the outflow boundary is shown. The 
eigensolution corresponds to 9 y = j and the eigenvalue /x = 0.8464 — 0.1382*. 


0.2 



-o® 1 1 1 1 

1800 1850 1900 1950 2000 2050 

i 

X 

Figure 11. Eigensolution on a grid with N x = 2048; a = 45°; 0 y = ^ ; and the corresponding eigenvalue is n = 0.8464 — 0.1382i. 

Even for combinations of grids and solutions with fast asymptotic convergence, many slow DC iterations may 
be encountered before the asymptotic regime is attained. Algorithmic enhancements are required to accelerate slow 
iterations preceding the asymptotic convergence and to improve asymptotic convergence, if necessary. Multigrid 
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addresses both these issues. Convergence of standard full-coarsening multigrid cycles for second-order convection 
discretizations on high-density grids is limited by the factor 0.75. 19 For the diagonal flow alignment the scheme 
becomes third-order accurate and the limiting factor is even more severe, 0.875. However, these rates are significantly 
better than slow-iteration DC rates. 

D. Defect-correction iterations for small-perturbation flows 



0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 

Mach 


Figure 12. Asymptotic convergence of DC iterations computed with the SD analysis. 



Figure 13. Asymptotic convergence of DC iterations computed with the SD analysis. 


In this section, DC iterations are applied for SP flows away from boundaries and singularities. The asymptotic 
rates of DC iterations are computed with the SD analysis. The angles of attack are varying as 0 < a < 45° and 
Mach number is varying between (almost) zero and fully supersonic, 0.01 < M < 1.81. Fig. 12 shows levels of the 
asymptotic rate on a 128 2 grid. The grid is not a high-density grid and the rates do not necessarily show the maximum 
values approached in grid refinement, but the distribution is representative for the medium-density grids. It shows that 
the slowest convergence is expected at low and near-sonic Mach numbers. 
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Fig. 13 shows the variation of asymptotic convergence rates versus Mach number on grids of progressively high 
density. The maximum rate over the range of angles of attack 0 < a < 45° is shown. Slowdown at low and sonic 
Mach numbers is observed on all grids. Similar to the convection convergence pattern shown in Fig. 10, the rates slow 
down for all Mach numbers on finer grids. 

Actual computations performed on the inflow/outflow domain shown in Fig 1 indicate similar trends. Fig. 14 
shows the “asymptotic” rate, namely, the last rate exhibited before achieving the machine-zero error, and the maximum 
convergence rate observed over the course of iterations. The rates shown in Fig. 14 are somewhat different from the 
rates predicted by the SD analysis because the error is sometimes reduced to the machine-zero level before the actual 
asymptotic convergence is achieved. As expected, the maximum rate is closer to one than the asymptotic rate. Both 
maximum and asymptotic rates peak at M ss 0 and M ~ 1. 



Figure 14. Convergence rates of DC iterations observed in actual SP flow computations. 
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